Toward a Consistent Prediction of Defect Chemistry in CeO2

Polarizable shell-model potentials are widely used for atomic-scale modeling of charged defects in solids using the Mott–Littleton approach and hybrid Quantum Mechanical/Molecular Mechanical (QM/MM) embedded-cluster techniques. However, at the pure MM level of theory, the calculated defect energetics may not satisfy the requirement of quantitative predictions and are limited to only certain charged states. Here, we proposed a novel interatomic potential development scheme that unifies the predictions of all relevant charged defects in CeO2 based on the Mott–Littleton approach and QM/MM electronic-structure calculations. The predicted formation energies of oxygen vacancies accompanied by different excess electron localization patterns at the MM level of theory reach the accuracy of density functional theory (DFT) calculations using hybrid functionals. The new potential also accurately reproduces a wide range of physical properties of CeO2, showing excellent agreement with experimental and other computational studies. These findings provide opportunities for accurate large-scale modeling of the partial reduction and nonstoichiometry in CeO2, as well as a prototype for developing robust interatomic potentials for other defective crystals.


INTRODUCTION
Ceria (CeO 2 ) is a technologically important rare-earth oxide with broad applications in several areas, including heterogeneous catalysis and solid oxide fuel cells (SOFCs). 1−3 It is versatile in catalytic applications because of its unique defect properties. Apart from its critical role in automobile three-way catalysts, ceria also holds great promise in novel catalytic processes such as selective methane oxidation to methanol, CO 2 conversion reactions, and hydrogen production through water splitting. 4−6 CeO 2 has outstanding oxygen storage capacity, which arises from its variable stoichiometry enabling the release or uptake of oxygen to pass through redox cycles. 7 As a consequence, reduced ceria (CeO 2−x ) is widely used as a supporting material for single-atom/nanocluster catalysts with exceptional catalyst-stabilizing and oxygen-releasing capabilities. 8−10 In addition, the high ionic conductivity of CeO 2 is necessary for its use as an electrolyte in SOFCs and relies to a large extent on the rapid oxygen-vacancy formation and migration processes. 11,12 Remarkably, CeO 2 remains very stable in the cubic fluorite structure over a wide range of temperatures. The earth-abundant nature of cerium 13 further makes it promising for large-scale industrial applications.
Oxygen vacancies are known to have low formation energies in ceria. 7 Accompanied by the removal of an oxygen ion to form a doubly ionized vacancy V O •• , two excess electrons could localize on two cation sites and form Ce Ce ′ small polarons 14,15 or be compensated by an oxygen interstitial O i ″ to form an anion−Frenkel pair. 16 Here, we employ the Kroger−Vink notation 17 for a point defect M S C , in which M is the defect species, S is the lattice site that the defect occupies, and C is the charge of the defect relative to the original site (dots, primes, and crosses represent positive, negative, and neutral charges, respectively). For the former process, Paier et al. 7 proposed that the best estimate of the charge-neutral oxygen vacancy formation energy in bulk ceria could be 4.2 ± 0.3 eV (relative to (1/2)E Od 2 ) in O-rich conditions deduced from conductivity and calorimetric titration experiments, close to the heat of reduction from CeO 2 to Ce 2 O 3 of 4.0 eV. 18,19 Surface vacancies have much lower (1−3 eV) formation energies than those in the bulk, depending on the facets and concentrations. 20 −24 In recent years, much effort has been devoted to studying the defect chemistry in CeO 2 by plane-wave density functional theory (DFT) techniques using periodic supercell approaches. 7,25−31 CeO 2 is a strongly correlated oxide. The description of the highly localized 4f states becomes problem-atic in DFT calculations when Ce 4+ is reduced to Ce 3+ . Standard Kohn−Sham DFT calculations using semilocal generalized gradient approximation (GGA) functionals such as Perdew−Burke−Ernzerhof (PBE) predict an unrealistic delocalized spin density for the occupied f states due to the self-interaction error (SIE). 25 While the Hubbard U correction scheme 32 (DFT+U, U Ce4f = 4.5−6.0 eV) can solve the electron localization problem, 7 the calculated formation energy of the oxygen vacancy (from 2.84 to 3.27 eV 26−29 ) is underestimated by 1.0−1.5 eV compared with experiment (4.2 ± 0.3 eV 7 ). Moreover, the choice of the U value is not straightforward, as it introduces uncertainty in the description of bulk properties, defect formation, and surface reactivity of ceria. 31,33 It has been reported that with the increase in U the lattice constant of CeO 2 becomes too large with increasingly underestimated oxygen vacancy formation energy. 31 Furthermore, Branda et al. 34 reported a variation of the oxidation state of Au on CeO 2 (111) with the U parameter. Watson et al. 25,26 proposed employing the DFT+U correction also to the O 2p states, which further properly localizes holes in ceria trapped by native defects and impurities. However, this approach underestimates the formation energy of the oxygen vacancy even more: 2.22 eV based on U Ce/O = 5.0/5.5 eV compared with 2.60 eV using U Ce/O = 5.0/0 eV. Alternatively, hybrid functionals that include a fraction of the nonlocal Fock exchange, although computationally more demanding, yield results in much closer agreement with experiment. Several hybrid functionals have been used to study the bulk and defect properties of CeO 2 based on periodic models. 30, 31,35−37 Notably, the oxygen vacancy formation energies predicted by hybrid functionals are much closer to experiment than the DFT+U predictions: 3.84 eV using PBE0 38 and 3.63−4.09 eV using HSE06. 7,30,31 Such improvements from hybrid functionals were also seen in other metal oxide systems with accurately predicted defect processes. 39 The configurations of oxygen vacancies and the associated Ce Ce ′ sites in reduced CeO 2 have been a topic of considerable debate. Early theoretical studies assumed that the two excess electrons are trapped at the nearest neighbor (NN) sites of the vacancy (i.e., NN-NN). 20, 25,40 Later, DFT+U calculations performed by Wang et al. 41 favored the localization of both electrons at the next-nearest neighbor (NNN) sites (i.e., NNN-NNN). In contrast, Murgida et al. 29 obtained three configurations (one NN-NNN and two NNN-NNN configurations) sharing the same lowest formation energies. Despite the minor differences, these predictions are generally consistent. Early molecular-mechanical (MM) Mott−Littleton (M-L) calculations on doped ceria also suggested that large trivalent cations prefer to locate at the NNN sites over the NN sites. 42 Scanning-tunneling microscopy (STM) imaging over the CeO 2 (111) surface conducted by Jerratsch et al. 43 captured several different configurations, indicating that at least one Ce Ce ′ is not adjacent to the vacancy site. The polarization energy of V O •• is large due to the high dielectric constants of CeO 2 , which can stabilize isolated charged species with respect to the neutral lattice, including both V O •• and Ce Ce ′ . Besides, the formation of Ce Ce ′ requires an expansion of the surrounding lattice, which opposes the location next to the oxygen vacancy site where the neighboring Ce−O bonds shrink. As a result, the binding energy of V O •• with the NN Ce Ce ′ polaron is very low (ca. 0.1 eV as reported by Sun et al. 30 ), indicating that Ce Ce ′ polarons are not tightly bound to oxygen vacancies. Experimentally, early conductivity measurements by Tuller and Nowick 44 on single-crystal CeO 2−x samples also revealed that V O •• is the predominant charge state at a small x in CeO 2−x (x < 10 −3 ), with a transition to singly ionized V O • at a larger x in more reducing environments. In addition, a recent cathodoluminescence spectroscopy study by Thajudheen et al. 45 demonstrated that the relative concentration of different charge states of oxygen vacancies depends on the oxygen partial pressure. Luo et al. 16 observed using neutron scattering that the anion−Frenkel pair is dominant in the bulk, while Ce Ce ′ polarons tend to aggregate at the surface and form an ordered reduced phase from nanorod samples. Such complexity in the charge compensation mechanism of oxygen vacancies could have a substantial impact on the ionic conductivity in SOFCs 12 and surface reactivity when CeO 2 serves as a catalyst 46 or support. 47 To resolve these uncertainties, reliable theoretical techniques are required to study the environmentally dependent formation mechanisms of all the relevant charge states of intrinsic defects in CeO 2 .
While being well-suited to modeling bulk materials and delocalized states, periodic boundary conditions have inherent limitations in modeling localized states such as isolated defects, polarons, and adsorbed molecules due to spurious image− image interactions. The long-range Coulomb interactions between periodic images originating from the net dipole moment of localized states are non-negligible. 48 As a result, periodic DFT calculations on defective solids usually require a large supercell to mitigate these spurious periodic interactions, 7,29 and additional correction schemes are usually necessary to achieve a higher level of accuracy. 48−50 Also, for some highly charged defects, the finite-size supercells could still be insufficient to fully consider the long-range atomic perturbations. In particular, the formation of highly charged defects in CeO 2 , such as Ce i •••• and V Ce ⁗, could provide a major perturbation of the surrounding atomic structure. These errors could be magnified when defect clusters are considered. Additionally, it is computationally very demanding with hybrid functionals to employ a sufficiently large supercell to minimize these errors.
A minor change in the formation energy of a point defect could result in several orders of magnitude variations in its calculated concentration, which could strongly affect the predicted electronic and optical properties. 39 Therefore, minimizing errors in the modeling method is necessary to ensure the accuracy of defect predictions. Embedded-cluster techniques that naturally avoid periodic image interactions and further consider long-range polarization effects are an effective approach to modeling localized defects in solids at the dilute limit. 51−56 The Molecular Mechanical, Mott−Littleton (MM M-L) approach 57−59 and hybrid quantum mechanical/ molecular mechanical (QM/MM) embedded-cluster techniques 55,60,61 are widely used for modeling defect processes in solids. Both methods are based on the embedded-cluster framework: the QM/MM model includes a QM core, an interface, and part of MM atoms in the active region, while the M-L model only has one active MM region, which is allowed to relax fully in response to the formation of charged defects; the outer fixed part reproduces the infinite bulk environment. 59 The M-L approach treats the outer part within a harmonic approximation to reproduce the bulk crystal field and linear dielectric response, while our QM/MM model implements a simplified continuous dielectric medium approximation where the response to the defect charge is calculated a posteriori. 55 Furthermore, the shell-model 62 interatomic potential (IP) is typically used in both techniques, in which the contributions of electronic long-range polarization effects due to the defect formation to the total energy can be computed through shell relaxations. With the rapid development of high-performance computing platforms, embedded-cluster calculations routinely employ over one thousand or more atoms in the active region (including QM and MM atoms) to consider explicitly the atomic displacements caused by defects, ensuring accuracy in predicting the formation energies and spectroscopic features. 54,63 M-L calculations of defects in solids are performed at the pure MM level, which relies heavily on the accuracy of the IP. A reliable shell-model IP is a fundamental prerequisite for quantitively predicting the defect structures and formation energies. In contrast, embedded-cluster results commonly use DFT calculations employing hybrid functionals (although higher level theory may be used 64−66 ) and with the IPpredicted dielectric response. Hence, an intrinsically consistent prediction of defect structures and formation energies by both approaches using the same potential is the ultimate goal that demonstrates the robustness of an IP. In previous work, such consistent predictions were achieved for MgO 67,68 and ZnO 51,69 but were limited to a few charge states of point defects without holes or electrons.
To date, a number of shell-model potentials have been developed for CeO 2 , which are collected in Table S1. 42,70−76 IP-based computational studies of defect chemistry in ceria were first performed by Butler, Catlow, Cormack et al. 70,77,78 using the M-L approach; these studies calculated vacancy migration energies in good agreement with experiment and demonstrated the crucial role of the ionic radius of dopants in determining the solution energy and the magnitude of dopant−vacancy interactions. Periodic models were also exploited to investigate the oxygen vacancy formation on surfaces and its role in CO oxidation. 71,79,80 Despite these advances, no single IP can offer a consistently accurate description of all the physical and chemical properties of CeO 2 , as shown in Figure 1. It was possible to reproduce accurately the main bulk properties of ceria, including phonon frequencies and the oxygen vacancy migration barrier based on a more complex potential model (IP9 in Table S1), 76 but the defect formation and surface energies proved to be significantly overestimated. One of the critical reasons for the less satisfying performance of previous potentials is the lack of reliable reference data in the IP development that exacerbates the errors in empirical parametrization.
In this work, we propose a novel methodology for IP development assisted by the QM/MM approach, which allows us to obtain accurate reference data for ionic polarizabilities, structure and formation energies of polarons, and intrinsic defects in dielectric materials. A new shell-model IP has been developed for CeO 2 that reproduces accurately the experimental structure, elastic and dielectric constants, defect and surface properties, and phonon dispersion. This approach also provides references for parametrizing localized holes and polarons, endowing powerful capacities for modeling the localization of charge carriers accompanied by defect formation. Based on the new IP, the calculated structure and Figure 1. Performance of the previous (IP1-IP9) and our newly developed (IP10a and IP10b) shell-model potentials for modeling CeO 2 , compared with experimental and first-principles results. IP10a is developed based on the pure Buckingham potential in short-range interactions, while IP10b is the refined final version with improved performance, using a more complex potential form. The calculated values of each observable based on the corresponding potential are shown in the blocks. The blocks are colored according to the relative errors of the predicted properties compared with these experimental or theoretical references. Lattice constant a 0 (extrapolated to 0 K), 33,81 elastic constants C 11 , C 12 , and C 44 , 82 phonon frequencies F 1u (TO), F 2g , and F 1u (LO), 82 static and high-frequency dielectric constants ε 0 and ε ∞ , 83 bulk modulus B 0 , 84  formation energies of charged defects by the M-L and QM/ MM approaches using hybrid functionals achieve a very high level of consistency, which is greater than that achieved in previous work. Our potential is parametrized entirely using the shell-model and pairwise potentials without complex manybody interactions, ensuring excellent computational efficiencies in studying complex defective systems. Moreover, the newly proposed strategies for developing reliable shell-model potentials for accurate defect predictions assisted by QM/ MM calculations could be extended to other systems.

METHODOLOGY
This work combines various computational techniques including plane-wave DFT calculations, IP-based lattice and M-L defect calculations, and the hybrid QM/MM embedded-cluster approach.

Plane-Wave DFT Calculations.
Plane-wave DFT calculations were performed at several levels of theory using the Vienna Ab-initio Simulation Package (VASP) 86 to set theoretical references for developing new IPs. Full computational details for calculating the bulk and surface properties of ceria are given in Section S1.1 of the Supporting Information.

Interatomic-Potential-Based Calculations.
The development of IPs and IP-based calculations are performed based on the Born model of ionic solids 87 using the General Utility Lattice Package (GULP) code. 88,89 In traditional shell-model potentials, pairwise interactions are described by the Buckingham potential: where r ij is the distance between two interacting ions and A, ρ, and C 6 are the parameters. Ionic polarizability is treated by the shell model, in which the massless shell is connected to the atomic core by a harmonic spring with a spring constant k 2 : The sum of the core and shell charges on an ion equals its formal charge. The ionic polarizability α in vacuum in the shell model is given by where Y is the shell charge. The electrostatic Coulomb interaction is calculated as where k e is the dimensional Coulomb constant (14. where A, ρ, C 6 , E, D, a, r 0 , and c 0 are the parameters of the potential. IP-based defect calculations were undertaken based on the M-L approach. 57,58 A point defect or defect cluster is embedded in the central part of the model (region I), where ionic relaxations and polarization are allowed explicitly. The surrounding regions (region IIa and IIb) that model the infinite solid are described within the linear response approximation. Region IIa acts as the interface, in which interactions with region I ions are calculated by explicit summation, but ionic displacements are calculated self-consistently using a harmonic approximation for the defect energy expanded around defect equilibrium configuration, as first proposed by Norgett 90 and implemented in the GULP code. 88 Region IIb extends to infinity to reproduce the periodic electrostatic environment of the solid, in which the long-range polarization energy of a charged defect outside the active region I is calculated using the macroscopic dielectric response of a perfect crystal (without a defect). The radii of region I and region IIa were set to 25 Å and 40 Å, respectively, which proved to be sufficiently large for considering the long-range polarization effects and provide well-converged defect energies (<0.1 eV).
Details of IP-based surface and vacancy migration calculations are described in Sections S1.2 and S1.3 of the Supporting Information.

QM/MM Calculations.
We employ the hybrid QM/MM embedded-cluster approach to model defect formation in bulk CeO 2 as implemented in the ChemShell 60,91 code. Our hybrid QM/MM embedded-cluster model is divided into five regions. In the QM region, electronic-structure calculations are performed at the hybrid DFT level using NWChem. 92 We have tested the performance of several pseudopotentials and basis sets to balance the computational accuracy and efficiency in large-scale QM/MM calculations of CeO 2 (with 100−200 QM atoms). The combination of the Def2-TZVP basis set 93 for O and a [4s4p2d3f ] basis set developed by Erba et al. 94,95 based on the Stuttgart−Dresden quasi-relativistic small-core (28 core electrons) effective core potential (ECP) 96 for Ce was determined as the best choice. The differences in the calculated formation energies of defects and ionization potential between the current setting and more complete basis sets are less than 0.1 eV. We utilized three hybrid functionals for defect calculations: hybrid GGA functional B97-2 97 and PBE0 98 and hybrid meta-GGA functional BB1K, 99 which include 21%, 25%, and 42% HF exact exchange, respectively. For consistency, we selected PBE0 results as the reference data in the IP development, together with those from VASP calculations using the same functional.
The MM calculations are performed using GULP. The MM part of the QM/MM structural model is divided into two regions: the MMactive region, which is allowed to relax during the calculations, and the MM-frozen region, which is fixed to reproduce the effect of the bulk environment. The outermost layer of the entire model includes point charges, which were fitted to eliminate the effects of surface termination and reproduce the Madelung potential of CeO 2 .
The interface region participates in both QM and MM interactions, serving as the buffer layer to minimize the mismatch of the QM and MM levels of theory. We used a specially designed local pseudopotential on the cationic sites in the form of a linear combination of three Gaussian functions. 56 The fundamental idea and detailed procedures of the QM/MM interface treatment are shown in Section S1.4 of the Supporting Information. The parameters were trained with a least mean square procedure for residual gradients on the ions in the active (QM, interface, and MM-active) region and the scatter of deep core levels (in this study, O 1s ) in the Kohn−Sham spectrum. The fitted pseudopotential for the interface cations has the form of QM/MM defect calculations were performed with the pythonbased version of ChemShell (Py-ChemShell) 91 on a large embeddedcluster model with ∼10 000 atoms in total and an O-centered 197atom QM cluster focusing on the formation of oxygen vacancies. The accurately calculated energies of polaron formation at the NN and NNN sites have been used for validation and refinement of the interatomic potentials that we report. A Ce-centered 111 QM-atom model was used for other types of defects, which yields convergence of defect formation energies to ca. 0.1 eV.

Calculations of Defect Formation Energies.
Lattice energy is the energy of the ionic compound with respect to constituent ions in the gas phase. The lattice energy ΔU L can be calculated according to the Born−Haber cycle 100 from experimental thermodynamical data, which (for CeO 2 ) is given by where ΔH subl (Ce) is the sublimation enthalpy of Ce (4.380 eV 101 103 obtained an average value of 9.41 eV from several oxides. Previous work has shown that ΔU L is a critical quantity that strongly affects the accuracy of M-L defect calculations. 104 However, fitting the potential to reproduce the "experimental" ΔU L according to an arbitrary choice of A O 2 may be problematic and could lead to several electronvolt errors in the calculated formation energies. We propose a self-consistent approach to determine the value of A O 2 in CeO 2 (8.14 eV, and therefore ΔU L = −107.5 eV according to eq 7) based on QM/MM calculations of defect energies and plane-wave DFT calculations of surface energies, which will be discussed in the following sections. For the defect formation energy calculations using the Born−Haber cycle based on the M-L results, a consistent usage of the A O 2 and ΔU L predicted by the IP gives accurate results that are comparable with QM/MM calculations.
It is necessary to clarify the definition of some energy terms in defect calculations. In QM/MM calculations, the formation energy of a defect in the charge state of q is defined as where E[X q ] and E 0 are the calculated total QM/MM energies of the defect and perfect structures, n i is the number of species that have been added (n i > 0) or removed (n i < 0) from the system to form the defect, μ i is the chemical potential of the species i, and E F is the Fermi energy relative to the valence band maximum (VBM). In order to account for the long-range polarization effect for charged defects that extends to infinity, an a posteriori correction is applied using the Jost formula, 61,105 where R is the radius of the active region (including QM, interface, and MM-active regions) and Q is the net charge of the defect system. High frequency ε ∞ and static ε 0 dielectric constants are used for vertical and adiabatic processes, respectively. M-L calculations compute the defect energy E[X q ] taking the gas phase ions (O 2− (g) and Ce 4+ (g)) as the reference, which differs from the formation energy (with respect to O 2 (g) and Ce(s)) in conventional DFT calculations. The energy of O 2 (g) and Ce(s) is not defined in IP-based calculations. Instead, such energies can be in turn obtained from the Born−Haber cycle 100 Therefore, the energy of the added or subtracted species in defect formation can be obtained by  100 Defect energies calculated by the M-L approach implemented in GULP have already accounted for the long-range polarization effect, so extra corrections are not needed.
For example, when an oxygen vacancy V O q with a charge state q (q = 0, +1,+2) is formed in CeO 2 , the excess 2 − q electrons are not trapped at the vacancy site but will localize on the nearby Ce site to form Ce Ce ′ small polarons. This process can be formed thus: The formation energy determined using M-L calculations is then calculated as the corresponding reaction energy: The formation energy of a charged defect is also dependent on the growth conditions. In the O-rich/Ce-poor conditions, the upper limit is determined by the formation of O 2 molecules, in which μ O = 0 eV and μ Ce = ΔH f (CeO 2 ). In O-poor/Ce-rich conditions, the lower limit of μ O is determined not only by CeO 2 but also by other Ce n O m phases, where CeO 2 should remain more stable than any other reduced phases. Conventionally, we take Ce 2 O 3 as the limiting phase as employed in previous DFT calculations. 27,28,30 Hence, Therefore, the chemical potentials should satisfy In both M-L and QM/MM calculations, formation energies (per defect) of the unbound cation−Frenkel pair, anion−Frenkel pair, Schottky trio, and interstitial disorder defects were calculated by

Strategies for Developing Robust Shell-Model Potentials.
In this work, we report a framework for developing a robust shell-model IP for CeO 2 assisted by QM/MM calculations, which could also be extended to other ionic solids. First, a classical shell-model IP based on Buckingham potentials was parametrized to reproduce the bulk structure and dielectric constants. As discussed in the previous section, to make the optimal choice, we systematically reviewed the performance of several candidate potential parameter sets based on different values of the ρ parameter in modeling defects and surface properties. Next, the optimal potential, IP10a, was used in hybrid QM/MM embeddedcluster models to obtain accurate references of the structure and formation energies of intrinsic defects in CeO 2 , including oxygen vacancies, cerium vacancies, oxygen interstitials, cerium interstitials, hole polaron, and electron polaron. Then, these data were used in turn to optimize the Ce 4+ −O 2− interactions using a more complex form that combines the Buckingham, Lennard-Jones, and Morse potentials for correcting the derivatives and an offset constant for correcting the lattice energy, IP10b. Finally, the Ce 3+ −O 2− and Ce 4+ −O 1− potentials were fitted based on the Ce 4+ −O 2− potential according to QM/MM calculations of electron and hole polarons in terms of bond length and formation energy. Based on the newly developed IP10b, the calculated defect structures and formation energies at the MM level of theory using the M-L approach are highly consistent with QM/MM results using hybrid DFT functionals.
3.1.1. Construction of the Reference Data Set. CeO 2 is a widely studied material with well-documented experimental measurements of physical and chemical properties. However, due to the easily reducible nature, oxygen vacancies and electron polarons in CeO 2 could significantly affect the measurements, resulting in inconsistent reports of several properties. For example, a considerable concentration of Ce 3+ (12%−45%) was detected in many experimentally synthesized polycrystalline thin-film samples, which are very difficult to eliminate even with oxygen plasma treatment at high temperatures. 106,107 Such a high concentration of Ce 3+ in CeO 2 could originate not only from surface defects but also from amorphous Ce 2 O 3 formed at the grain boundaries of nanocrystals. 106 Hence, we focus on data obtained from single   Gupta and Singh reported an a 0 of 5.401 Å at 100 K, and no data are available close to 0 K. Based on several thermal expansion experiments, we obtained an estimate for a 0 of 5.395 Å when extrapolated to 0 K. 81,111−113 This estimate also allows us to obtain an a 0 of 5.411 Å at 300 K based on the new IP using free energy calculations within the quasi-harmonic approximation (20 × 20 × 20 k-points on the 12-atom conventional cell of CeO 2 ), in excellent agreement with experiment.
The O 2p → Ce 4f band gap of CeO 2 is typically reported ranging from 3.0 to 3.5 eV 3,106,114−116 and is known to decrease with the increase in the concentrations of oxygen vacancies and Ce 3+ . 3,106 However, Pelli Cresi et al. 117 recently proposed a revised value of 4 eV for the optical band gap based on steady-state and ultrafast transient absorption spectra measured on stoichiometric samples. They argued that the long absorption tail from 3 to 4 eV in the Tauc plot should be identified as the Urbach tail, further confirmed by the photobleaching in the transient absorbance measurement, which corresponds to the occupied Ce 4f localized state due to the rapid small polaron formation. Such an Urbach tail is also seen in nanoparticle samples with a larger band gap of 4.2−4.3 eV due to the quantum confinement effect. 118 This conclusion is also supported by previous high-resolution electron energy loss spectroscopy (EELS) measurements on stoichiometric 119 and fully oxidized single-crystal 120 samples, which located the unoccupied Ce 4f state 4−4.4 eV higher than the O 2p state. Hence, we consider that a band gap of 4 eV measured by Pelli Cresi et al. 117 could be the best estimate for CeO 2 and is further employed in this work. We will compare with experimental data for the elastic constants C 11 , C 12 , and C 44; 82 zone-center phonon frequencies F 1u (TO), F 2g , and F 1u (LO); 82 static and high-frequency dielectric constants, ε 0 and ε ∞ ; 83 bulk modulus B 0 ; 84 and migration barrier of oxygen vacancy E Vd O migration . 85 Plane-wave DFT calculations using VASP at different levels of theory were performed alongside our semiclassical GULP simulations. The calculated lattice constants, band gaps, Born effective charges, and surface energies of CeO 2 are listed in Table 1, where we find that the hybrid functionals in our assessment set show a superior agreement with experiment. Previous work also showed that in general hybrid functionals perform much better than GGAs in describing the formation enthalpy, heat of reduction, and defect formation energy of CeO 2 . 7,30,36 Because the PBE0 functional is implemented in both the plane-wave DFT code (VASP) and atomic-centered basis set code (NWChem), employed in our hybrid QM/MM calculations, we consistently use PBE0 results as the primary reference in the IP development. The Born effective charges and surface energies from VASP calculations, as well as the inlattice ionic polarizabilities, defect structures, and energies from QM/MM results, were used in refining the IP. Our new IP can easily be reparameterized to reproduce the predictions from any other functional based on the same protocol.
3.1.2. Strategies for Improving the Accuracy of Shell-Model Potentials. A shell-model potential can be divided into two parts: the shell model (Y and k 2 ) and other parameters including short-range repulsion and dispersion interactions. The shell-model defines the ionic polarizabilities that significantly affect the dielectric properties of materials, which play a key role in defect formation. Previous parametrizations of interatomic potentials for ceria have in common a somewhat unsatisfactory approach to choosing how the lattice polarizability is distributed between individual ions, which may not correctly reflect the physical nature of the target material. Polarizabilities have often been determined by fitting spring constants and shell charges to experimental data without consideration of the relative polarizabilities of cations and anions; it is possible to obtain a number of candidate potential sets with different shell-model parameters but with similar performances in modeling the bulk properties. An alternative to purely empirical fitting was developed by Lewis and Catlow 74 considering in-lattice polarizabilities based on gasphase Pauling's polarizabilities of cations, which are known to be less affected by the crystal environment. By contrast, the anion polarizability varies strongly with the oxide structure and chemical composition. 74 The in-crystal ionic polarizabilities can be calculated from hybrid QM/MM embedded-cluster models (cf. ref 127). Here, we have employed ChemShell with NWChem as a QM driver to construct suitable embeddedcluster models and calculated the in-lattice ionic polarizability. The detailed methodology is presented in Section S1.5 of the Supporting Information.
In Table 2, we collected the calculated in-lattice polarizabilities of Ce 4+ and O 2− ions using the Def2-TZVP 93 basis set. It can be seen that the gas-phase cation polarizability is very close to the frozen in-crystal polarizability, as has been pointed out by Lewis and Catlow. 74 The gas-phase calculations for O 2− are extremely basis-set-dependent, (as expected, since the O 2− ion is an unbound species in the gas phase), and increases from 6.293 Bohr 3 using the Def2-TZVP basis set to 50.871 Bohr 3 using Def2-QZVPPD 128 at the PBE0 level of theory. In the limit of a complete basis set, the ground-state solution would include an electron dissociated from the oxide ion to infinity and, therefore, an infinite polarizability. The hybrid DFT predictions are in good agreement with higherlevel coupled-cluster results. The intrinsic polarizability of an ion, α i , in the CeO 2 crystal environment does not include effects of charge transfer from or to surrounding ions, which plays a significant role in real materials, and the total polarizability of an ion is = + i ct (27) where α ct is the respective contribution from charge transfer. When fitting IP parameters, we kept the ratio of Ce 4+ and O 2− polarizabilities (α = Y 2 /k 2 ) constant according to the QM/MM calculated polarizability data, based on a conjecture that α ct is proportional to α i . Under this methodology, a unique set of shell-model parameters can be determined that models more Polarizability also determines the dispersion (London) interaction, the strength of which is usually described by the C 6 coefficient in the Buckingham interatomic potential. The similarly polarizable Ce 4+ and O 2− in CeO 2 suggest that the C 6 coefficients for Ce 4+ −O 2− , Ce 4+ −Ce 4+ , and O 2− −O 2− should all be considered, whereas most of the previous IPs did not introduce the C 6 coefficients for all these interactions. We have calculated C 6 according to the Slater−Kirkwood formula 129 using the participation numbers reported by Pyper et al. 130 and our calculated polarizability data. C 6 (XY) is given as where P X is the participation number of ion X (4.455 for O 2− and 7.901 for Ce 4+ from Pyper et al. 130 ).
In ionic solids, the long-range electrostatic interaction contributes greatly (over 90%) to the total lattice energy. Although the short-range repulsive interaction originating from the overlapping electron densities contributes much less, the pairwise parameters in an IP largely control the bond lengths and structures and, hence, energies and physical properties. Among all the short-range interactions, the cation−anion interaction has the greatest effects due to its directly bonded nature. The cation−cation and anion−anion interactions are much weaker because of the longer interatomic distances. We note that parameter ρ in the previous Ce 4+ −O 2− Buckingham potentials varied from 0.335 Å to 0.429 Å. A small ρ with a large A or a large ρ with a small A may show only minor differences in the calculated bulk properties. However, we found that the choice of ρ is of vital importance in predicting defect and surface properties, which are typically not included in the list of fitted observables. Hence, the "best" IP derived by the fitting procedure may not be appropriate for modeling defects and surfaces, as can be seen in the predicted incorrect relative stabilities of surfaces and defect pairs by several IPs.
To explore the effects of ρ, we parametrized several IPs based on different fixed values of ρ within the traditional Buckingham potential framework (shown in Table S2). These potentials share the same Ce 4+ −Ce 4+ and O 2− −O 2− parameters, while other parameters are fitted to minimize the sum of squares of selected observables. We used the A and ρ parameters of the O 2− −O 2− Buckingham potential from Lewis and Catlow 74 since they have remarkable transferability in many ionic oxide systems and can be combined with other potentials to study dopants and solid solutions. We also put the highest fitting weights on the lattice constant (ca. 5.395 Å at 0 K) and static and high-frequency dielectric constants due to their critical importance in determining the accuracy of QM/ MM calculations. Therefore, all the IPs we derived give accurate reproductions of the lattice constant and dielectric constants. Figure S1 shows how the ρ parameter affects the predicted properties of CeO 2 . The dark gray horizontal lines represent the reference data from experiment or DFT calculations at the PBE0 level of theory. If we only consider the bulk properties that are normally used as fitting observables, a larger value of ρ, such as ρ = 0.37 Å, could be the result obtained by the leastsquares procedure. However, if we further consider the defect and surface properties, those IPs with a larger value of ρ produce inaccurate results for ceria. By systematically considering bulk, defect, and surface properties (discussed in Section S1.6 of the Supporting Information), we determined the IP with ρ = 0.349 Å as the best Buckingham potential for CeO 2 , which was further used in QM/MM defect calculations. We stress that the QM/MM calculations based on IP10a are reliable and consistent with those based on IP10b since the lattice constants, dielectric constants, and ionic polarizabilities predicted by the two potentials are identical. The full set of parameters is provided in Table S2, and the performance in reproducing the properties of ceria has also been shown in Figure 1, denoted as "IP10a". Within the current model, it is still impossible to reproduce simultaneously all targeted properties because the cation−anion short-range interaction Chemistry of Materials pubs.acs.org/cm Article cannot be accurately described by a single Buckingham potential. However, we have shown here that the shell-model potential with the classical Buckingham form can qualitatively and semiquantitatively reproduce the correct atomic structure and properties of CeO 2 based on a careful selection of some critical parameters.

Refining the Interatomic Potential According to QM/MM Calculations.
Short-range interactions in solids sometimes cannot be exactly described by a single potential adopting a simple analytical form. The Lennard-Jones potential is known to be considerably too repulsive at short bond distances. 131 Although the Buckingham (and its constituent Born−Mayer) potential is a good physically justified approximation for the short-range repulsion, 132 empirical fitting procedures could cause unwanted problems, as illustrated by the critical choice of ρ described above. Combining different forms, nevertheless, could increase the parametric space and overcome the weakness of each component to improve the accuracy. 104,133 Derivatives of an IP determine the structure and elastic, dielectric, and phonon properties of the modeling system. The absolute potential value determines the lattice energy, which further affects the calculated energies of defects and surfaces. To obtain a robust IP, both derivatives and absolute values of the potential should be accurately determined.
We define a new analytic form for the cation−anion pairwise potential, as shown in eq 5. Parameters of the refined potential (IP10b) for CeO 2 based on this form are shown in Table 3. The sole Buckingham potential representing the short-range interaction is insufficient for a consistently accurate description of elastic, dielectric, and phonon properties. This weakness can be attributed to inaccurate potential derivatives near the first and second neighbor bond distances. Hence, a Morse-like function is superimposed on the Buckingham potential to correct these derivatives and reproduce these bulk properties accurately. The offset constant c 0 shifts down the potential curve by a constant to correct the lattice energy. Weak Lennard-Jones potentials were used to add repulsion to the interactions at very short bond distances to overcome the Coulomb catastrophe that could occur in molecular dynamics or Monte Carlo simulations. 133 The overall potential is truncated at 5.278 Å where the potential energy drops to 0 eV, between the second and third Ce−O neighbor bond lengths. The resulting potential curve is shown in Figure S2, compared with previous potentials in the literature. The shell model parameters were slightly optimized for better overall performance, while the ratio of polarizabilities of Ce 4+ and O 2− remained consistent with IP10a.
We here also emphasize the critical role of the lattice energy predicted by a shell-model potential in determining the accuracy of defect and surface calculations. The lattice energy ΔU L sets an essential basis for every calculation regarding energy differences such as surface energy, defect formation energy, and migration barrier. Lattice energy is not a direct experimental observable but can be evaluated theoretically through the Born−Haber cycle based on the experimental formation enthalpy, as shown in eq 7. However, the exact value of ΔU L requires the knowledge of the oxygen second electron affinity A O 2 , which is a crystal-structure-dependent variable 134 and not known from experiment. We proposed that the required ΔU L and A O 2 in the IP model can be determined according to DFT-calculated surface and defect energies. As shown in Figure 2, several variants of potentials are derived based on different ΔU L values by only varying the c 0 term while keeping all other parameters constant. As a result, these potentials share the same shape (derivatives) but differ in the absolute potential energy, thus predicting the same equilibrium-state structure but different absolute energies for defects and surfaces. The predicted defect and surface energies change considerably despite the minor differences in the predicted ΔU L . Periodic DFT calculations of surface energies and QM/ MM calculations of defect energies at the PBE0 level were employed as references to obtain the required lattice energy that gives the most accurate prediction. In IP-based surface calculations, the lattice energy mainly determines the required energy to generate the unrelaxed surfaces rather than the  (111) surface has only 0.01 (0.05) J m −2 relaxation energy as predicted by periodic DFT (IP-based) calculations, which could be the best reference for the target lattice energy. We found that both defect and surface calculations support that ΔU L = −107.5 eV, corresponding to an A O 2 of 8.14 eV, which is employed as the final version of our CeO 2 potential (IP10b). We note that the new potential still has certain errors on the calculated surface energies compared with periodic DFT results. To our knowledge, the only experimental measurement of surface energies of CeO 2 was based on nanoparticle samples (which should be dominated by the most stable (111) surface) that reported an average value of 1.20 ± 0.2 J m −2 , 125,126 slightly higher than DFT predictions. The overestimation mainly originates from the cleavage energy before structural relaxation, which is also seen in other IP models. We have checked that the optimized surface structures are consistent with PBE0 predictions with an average error of 0.012 Å for the predicted surface and subsurface Ce−O bond lengths, suggesting a reasonably good description by the IP model.

Parameterization for Hole and Electron Polarons in CeO 2 . Finally, we developed additional parameters for modeling localized electrons (Ce Ce ′ ) and holes (O O
• ) in CeO 2 . The electron small polaron has been experimentally observed by STM imaging and confirmed by DFT calculations on the (111) surface. 12,30,40,43 Holes also stabilize as small polarons in CeO 2 , 25 which is a common feature of many oxide materials. 135 The formation of electrons and holes could be modeled through the M-L approach as described by Freeman and Catlow in SnO 2 . 100 Conventionally, electrons/holes were modeled by adding/subtracting a unit charge from the original shell charge and keeping the same core charge. However, using the same short-range potential to describe interactions between ions in different oxidation states is less accurate and usually results in several electronvolt errors in the predicted ionization potential, electron affinity, and band gap compared with experiment or DFT. 100,104 Instead, we propose a correction scheme based on the new O 2− −Ce 4+ potentials but where we target the bond lengths and formation energies of small polarons obtained from our QM/MM calculations.
In both M-L and QM/MM calculations, vertical processes were modeled through electronic (shell) relaxation, while adiabatic processes required explicit structural optimization of the active regions. We performed QM/MM calculations to obtain the vertical and adiabatic ionization potentials and electron affinities. In the adiabatic processes, electrons and holes are self-trapped, forming small polarons in CeO 2 . Before fitting the potentials for modeling the polarons, the shell charges were first corrected in line with the formal charges of Ce 3+ and O 1− . Then, the coefficients in the Ce 4+ −O 2− potentials (A and D) are multiplied by a factor to correct the derivatives of the potential and reproduce the polaron structures. This factor was determined by targeting the average Ce 3+ −O 2− /O 1− −Ce 4+ first-neighbor bond length around the ionized species calculated by the M-L approach to correspond to PBE0 QM/MM results. Finally, the potential is offset by a c 0 constant to reproduce the polaron formation energy calculated by QM/MM. A similar procedure was used to parametrize the Ce 4+ −O 1− potential for describing localized holes, while the vertical ionization potential was used as the reference for correcting c 0 , ensuring the alignment of the positions of VBM in both techniques for accurate formation energy calculations of charged defects. Full details of the additional parameters have been given in Table 3.
In summary, we have provided some new strategies to develop shell-model potentials with improved performance in modeling ionic materials in Section 3.1. Implementation of inlattice ionic polarizabilities, lattice energy, defect energies, and surface energies in the potential fitting improves the description of various physical properties, while consideration of QM/MM results of charge carriers further makes it possible to model various charged defects at the MM level of theory. Based on this approach, the formation of native defects in different charge states can be reasonably described through the M-L approach, assuming that the charge carrier is well localized and trapped by defects, as is the case in CeO 2 . Such an approximation gives reasonable accuracy with results comparable with DFT calculations using hybrid functionals, as will be shown in Section 3.2.

Performance of the New Potential in Modeling
Defects in CeO 2 . The structure and properties of CeO 2 predicted by the revised potential have been shown in Figure 1, denoted as IP10b. Compared with IP10a, the revised version significantly improves the performance in many aspects: the predicted bulk modulus, elastic constants, and phonon frequencies are in excellent agreement with experiment, and the calculated defect energies are also more consistent with QM/MM predictions. Based on the new potential, we also calculated the phonon dispersion of CeO 2 at 0 K ( Figure 3). Experimental results are presented in markers for comparison, derived from reflectivity by inelastic neutron scattering by Clausen et al. 136 and Marabelli and Wachter 137 and Raman scattering by Nakajima et al. 82 and Kourouklis et al. 138 Our potential shows excellent agreement with the experimental measurements in terms of the phonon modes at the high symmetry points and acoustic mode dispersion, ensuring good accuracy for future study of thermodynamic properties of defective CeO 2 using free-energy calculations.
A detailed comparison of QM/MM and M-L calculations of the vertical and adiabatic formation energies of holes and electrons in CeO 2 is given in Table 4. Full computational details are given in Section S1.7 of the Supporting Information. For comparison, if we directly employ the Ce 4+ −O 2− shortrange potential to model the localized electrons and holes without any correction (listed in the first column), as originally proposed by Freeman and Catlow, 100 the bond lengths around the polarons and the band gap of CeO 2 are significantly overestimated. Our corrected potential accurately reproduces the formation mechanisms of localized electrons and holes in CeO 2 and gives a more reasonable band gap. This strategy could be a general approach and extended to other systems with localized charge carriers. The self-trapping energy of the localized electron is predicted to be −0.92 eV, compared with −0.61, −0.69, and −1.39 eV by the QM/MM approach using B97-2, PBE0, and BB1K functionals, respectively. Previous periodic DFT calculations reported −0.30 eV and −0.54 eV by HSE06 and PBE+U. 30 The predicted band gaps by the embedded-cluster approaches are higher than the experimental values for two reasons. First, despite the advantages in predicting bulk properties and formation energies of defects, most hybrid functionals significantly overestimate the band gap of CeO 2 . 36,37 CeO 2 is a strongly correlated f electron oxide, and the noncorrelated HF method predicts the band gap as 14.5 eV. 37 As a result, the hybrid functional that includes a higher percentage of HF exchange usually provides a worse description of the band gap of CeO 2 . However, a certain percentage of HF exchange (usually 25%) is required to ensure the predictions of localized small polarons in CeO 2 , or both electrons and holes will become more delocalized even when trapped by defects, which is however contrary to experiment. 7,25,43 Second, the overestimation of the band gap is partly due to the poor description of delocalized states and quantum confinement effects by the hybrid QM/MM embedded-cluster approach. Hence, it can be seen that the predicted "band gap" through I vertical − A vertical by the QM/MM method using the PBE0 functional (4.59 eV) is slightly higher than plane-wave periodic DFT calculations using VASP (4.38 eV) and the experimental measurement (4 eV). To counter these problems, we typically use the experimental band gap and calculated vertical ionization potential to determine the position of conduction band minimum (CBM) in QM/MM calculations. 56    In M-L calculations, the first column shows the results using modified shell charges for Ce 3+ and O 1− but the same potential for Ce 4+ -O 2− interactions, while the second column shows the calculated results using IP10b presented in Table 3. b Present work. c Reference 30. . Among these types of defects, O 2− ions are more easily affected by defect formation than Ce 4+ ions (except Ce i •••• , which is less common in CeO 2 due to the relatively high formation energy). A detailed comparison of the structural relaxation due to defect formation is made in Table S3. Overall, M-L calculations show excellent agreement with previous periodic DFT and our QM/MM predictions of the defect structures.

Chemistry of Materials
We also studied the formation energies of different charge states of native defects using the new IP and the embeddedcluster approaches. We started from the various possible configurations of polarons trapped by oxygen vacancies in bulk ceria. Previous DFT calculations have predicted that the oxygen vacancy becomes less stable when the two polarons locate beyond the second-neighbor sites. 29,42 This conclusion is confirmed by our M-L calculations. The calculated formation •• −2Ce Ce ′ ] × formed in bulk ceria, the excess electrons will localize on two Ce sites (Ce Ce ′ ). As shown in Figure 5 and Table 5, there are one NN-NN, three NN-NNN, and five NNN-NNN symmetrically unique configurations on which the two Ce Ce ′ polarons can locate. 41 The calculated formation energies of all nine configurations Table 5. In general, our simulations confirm that the excess electrons have lower energies when localized on the NNN site, consistent with other DFT predictions. 29,41,139 The most stable configuration is the NN1-NNN2 with a formation energy of 3.564 eV, and the second stable NNN1-NNN4 configuration is only 0.014 eV higher in the formation energy. Our QM/MM calculations also favor the stabilization of the NN1-NNN2 configuration against the NNN1-NNN4 configuration, consistent with the M-L predictions. The tiny energetic deviations (<0.12 eV) among the formation energies of all the nine configurations suggest very close stabilities of different configurations, consistent with the previous periodic DFT+U study by Murgida et al.. 29 Because the adiabatic hopping barrier for Ce Ce ′ in CeO 2 is very low (0.15−0.2 eV 12,140,141 ), the positions of polarons could change rapidly at elevated temperatures. Previous DFT calculations predict lower formation energy at the PBE+U level (2.2−2.95 eV) 26,29−31 and slightly higher formation energy using HSE06 (3.63−4.09 eV 7,30,31 ). Brugnoli et al. 38 reported 3.84 eV for the formation energy of the NN1-NN2 configuration based on the periodic model using PBE0, compared with 4.04 eV from our QM/MM model using the same functional. Burow et al. 142 conducted the periodic electrostatic embedded-cluster method (PEECM) using PBE0 but obtained a much lower formation energy (3.0 eV) than our QM/MM predictions and periodic calculations from Brugnoli et al. 38 Our predictions at the PBE0 level of theory are also in line with the experimental heat of reduction of ceria, which is ca. 4.0 eV. 18, 19 The same approach was used to investigate other charged defects with localized electrons and holes. Table 6 summarizes the calculated formation energies of all charge states of point defects studied by the M-L and QM/MM methods using our new IP. Figure 6 shows a direct comparison of the calculated formation energies under different conditions using the M-L approach and QM/MM calculations. These two embeddedcluster approaches present a high level of consistency in predicting the formation energies, confirming the success of our newly proposed fitting strategy for localized charge carriers. For the charge-neutral point defects compensated by electrons and holes, our predictions are in good agreement with periodic calculations. For highly charged defects, the longrange polarization effects should greatly contribute to the calculated formation energies as evaluated using the Jost correction in eq 9, which is −0.46 eV, −1.84 eV, −4.14 eV, and −7.37 eV for the ±1, ±2, ±3, and ±4 charge states, respectively. These effects were typically not explicitly considered in previous supercell models. Recently, there has been increasing awareness of the application of long-range electrostatic corrections toward more quantitative predictions  Our M-L predictions of the formation energies of chargeneutral defect pairs and trios are also in good agreement with QM/MM predictions. When an anion−Frenkel pair is formed, a lattice oxygen ion occupies an octahedral interstitial site, leaving an oxygen vacancy. The formation energy is predicted to be 2.15 eV per defect (4.30 eV for the total defect reaction) by our potential model at infinite separation (2.23 eV, 2.24 eV, and 1.77 eV at the QM/MM B97-2, PBE0, and BB1K levels of theory, respectively, 2.07−2.08 eV by periodic PBE+U calculations, 27,28 and 1.83 eV by experiment 144 ). In our M-L calculations, we observed the recombination behavior of bound anion−Frenkel pairs in the nearest (2.34 Å) and nextnearest (4.47 Å) separations. The first stable anion−Frenkel pair has a separation distance of 5.88 Å with a total formation energy of 3.85 eV (viewed as a defect cluster), yielding a binding energy of 0.45 eV. The second stable pair has a slightly higher formation energy of 3.98 eV due to the larger separation (7.01 Å). Our M-L calculations predict a similar critical stabilization distance of the bound anion−Frenkel pair as reported by Huang et al. 28 16 This behavior is similar to other cubic fluorite crystals such as UO 2 and cubic ZrO 2 . 72,145 Oxygen vacancy formation is of crucial importance in bulk ceria, which can repeatedly absorb and release oxygen in catalytic reactions through the Mars−van Krevelen mechanism. 9,10 Furthermore, the V O •• migration barrier is very low (0.53 eV predicted by our potential, 0.5−0.6 eV by previous DFT and experimental studies 30,85,146−148 ), which could be further reduced by the surrounding Ce Ce ′ polarons on the surfaces, 12,140 ensuring good ionic conductivity in SOFC applications. All these features make this material versatile in energy and catalytic applications.
The energetic difference in predictions between the embedded-cluster approaches and periodic PBE+U results could mainly originate from the differences in DFT functional and long-range polarization effects. As illustrated before, the periodic PBE+U approach may not be accurate enough for the formation energy of highly charged defects, while our QM/ MM calculations that employ hybrid functionals and consider long-range polarization should be more reliable and target the defect formation at the dilute limit. Overall, IP10b accurately describes defect structures and formation mechanisms in ceria, which could be useful for investigating more complex defect clusters and partially reduced phases that are computationally too expensive for first-principles calculations. Moreover, our potential model for modeling holes and electron polarons is highly transferable. Using the same correction approach, one can reproduce the predictions of other hybrid functionals. 2 . The formation energy of a point defect depends on the position of the Fermi level (electronic chemical potential) and chemical potentials of exchanged species according to eq 8, which will be further affected by conditions, especially oxygen partial pressure P and temperature. 107 As summarized by Reuter and Scheffler,149 the relationship between the oxygen chemical potential μ O (T, P) and oxygen partial pressure P can be obtained by

Thermodynamic Defect Processes in Bulk CeO
where P 0 = 1 atm is typically defined as the zero state, eV, H is the enthalpy, and S is Chemistry of Materials pubs.acs.org/cm Article entropy. Using the experimental enthalpy and entropy data, 150 we calculated the value of μ O at various reaction conditions. As shown in Figure 7a, μ O decreases with increase in temperature or decrease in P. In the case of CeO 2 , an extreme reducing condition (e.g., P = 2.3 × 10 −10 atm at 1800 K) is required to reach the O-poor limit of μ O min = −3.94 eV at which Ce 2 O 3 becomes more stable than CeO 2 . In realistic systems, the formation energies of defects always lie between the two limits shown in Figure 6. The Fermi level of CeO 2 is in turn also affected by temperature and oxygen partial pressure. 107 Variation of the Fermi level influences further the dominant defect type, energetics of charge states, and concentrations of defects in solids. The self-consistent Fermi energy (shown in Figure 7b) was determined using the SC-FERMI code 151 based on the M-L defect formation energies and density of states of the pristine CeO 2 (obtained using VASP with PBE0 functional, where the unoccupied states were shifted downward to reproduce the experimental band gap of 4 eV). In the O-rich limit, the Fermi level lies ca. 2.5 eV above the VBM and favors the formation of doubly ionized V O

••
. These vacancies could be compensated by acceptor defects such as O i ″ and O i ′ or electrons that are not trapped at the vacancy sites. With the increase in temperature or decrease in oxygen partial pressure, the Fermi level rises toward the CBM, which could further support the stabilization of V O • and V O × with trapped electron polarons. An electrical conductivity study by Tuller and Nowick 44 on single-crystal CeO 2−x samples also concluded that V O •• is dominant at small x in CeO 2−x (x < 10 −3 ), with a transition to singly ionised V O • at greater x under reduced P, confirming our predictions. The thermodynamic transition level ε(q 1 /q 2 ) is defined as the energy of an adiabatic ionization process changing the charge state q 1 of the defect to another charge state q 2 . It can be obtained as the electron Fermi level where the formation energies of two charge states of a point defect are equal or using the following formula: 39 where E f (X q , E F = 0) is the formation energy of the defect X with the charge state q when the Fermi level is at the valence band maximum (VBM, E F = 0). The thermodynamic transition level indicates the relative stability of the two charge states. When E F < ε(q 1 /q 2 ), the charge state q 1 is stable and vice versa. The calculated thermodynamic transition levels of oxygen vacancy are shown as dashed lines in Figure 7b. In our M-L calculations, ε(+2/+1) = 3.48 eV and ε(+1/0) = 3.51 eV above the VBM. Sun et al. 30 have compared DFT PBE+U and HSE06 predictions using the supercell approach with an appropriate correction scheme on charged defects. They also obtained high transition levels for oxygen vacancies using the HSE06 functional (ε(+2/+1) = 3.1 eV and ε(+1/0) = 3.2 eV), but their PBE+U calculations yielded much lower values of ε(+2/+1) = 1.45 eV and ε(+1/0) = 1.65 eV. This difference could originate from the accuracy of the predicted band gap from the DFT functionals.
Our previous discussion focused on the equilibrium-state defect formation in bulk CeO 2 at the dilute limit, without considering the extrinsic doping and surface effects. Even trace dopants, or impurities can considerably affect the position of the Fermi level and defect formation. 26,70,152 Further, we should note that the formation energies of V O × on or near the surfaces are much lower than their counterparts in bulk, which could be rationalized by a decrease in the coordination of atoms exposed at the topmost surface layer. Previous PBE+U calculations reported an energetically favorable formation of V O × in the subsurface layer of CeO 2 (111) with a formation energy of ca. 1.8 eV compared with 2.85 eV in bulk. 24,153 The formation energies of V O × on less stable CeO 2 (110) and CeO 2 (100) surfaces are predicted to be even lower by Ganduglia-Pirovano et al. 24 and are 1.06 and 1.35 eV, respectively. Such low formation energies of V O × should result in nonnegligible shifts of transition levels toward the VBM and would favor the trapping of excess electrons near the surfaces, as experimentally captured by STM images. 40,43 This observation is, in fact, quite general. For example, recently, Wang and Yin 154 also found that iodine vacancies in CH 3 NH 3 PbI 3 formed in the bulk and on the surface have significant differences in the transition level and effects on photovoltaic performance. Returning to CeO 2 , reduction in practice always starts from the near-surface layers, resulting in the localization and segregation of excess electrons, which could lead to further reconstruction into alternative ordered phases such as CeO 1 . 6 8 , Ce 7 O 1 2 , Ce 1 1 O 2 0 , and Ce 3 O 5 . 16,110,155,156 Trace dopants or unintended impurities including divalent 157 and trivalent 158 metal ions and surface absorbates such as fluorine 159 and hydroxyls 160 could also become the charge-compensating species for V O •• and have a large impact on the Fermi level and defect concentrations. Our calculations show that, with the increase in temperature or decrease in oxygen partial pressure, the Fermi level rises in bulk CeO 2 , and singly ionized V O • and charge-neutral V O × become dominant before the formation of ordered reduced phases. A complete understanding of the partial reduction behavior of ceria requires further investigation using, e.g., free energy calculations, Monte Carlo, or molecular dynamics techniques 161−164 to consider explicitly entropic effects of defect formation and defect−defect interactions, which will be investigated in future research based on our new potential.

CONCLUSIONS
We have proposed a new strategy for deriving shell-model interatomic potentials based on hybrid QM/MM embeddedcluster calculations of ionic polarizabilities, defect structures, and formation energies. A new potential has been developed for CeO 2 , with a distinctive performance in predicting the structure, elastic, dielectric, defect, surface, phonon, and thermodynamic properties compared with experimental measurements and ab initio calculations. In particular, the calculated structures and formation energies of polarons, various charged states of native defects, and defect pairs in CeO 2 achieved a high level of consistency between the M-L calculations and QM/MM results. Benefiting from our new potential, one could provide more insights into the complex defect chemistry of CeO 2 . These developments will enable future work on defect formation at elevated temperatures using free energy calculations; partial reduction and nonstoichiometry based on molecular dynamics and Monte Carlo simulations; and the role of surface defects in catalytic reactions using QM/MM techniques. This novel strategy for developing robust shell-model potentials capable of predicting accurate defect properties could be further extended to other ionic systems.